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1.  Introduction 

Nonlinear,  parametrized  equations 
F(z,X)  =  0, 
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represent  models  of  equilibrium  problems  for  many  physical  systems.  If 
F=  Rn  -*  Rm  ,  n=m+p,  p  i  1.  is  continuously  differentiable  on  Rn.  then  the 
regular  solution  manifold 

n  =  {  x  e  Rn ;  F(x)  =  0  ,  rank  DF(x)  =  m  }  (1.2) 


is  a  p-dimensional,  differentiable  manifold  in  Rn  without  boundary.  We 
shall  assume  always  that  F  is  at  least  of  class  Cr,  r  i  2. 

The  standard  procedures  for  the  computational  analysis  of  such 
solution  manifolds  are  the  continuation  methods.  When  the  parameter 
dimension  p  exceeds  unity,  these  methods  require  a  restriction  to  some 
path  on  the  manifold  and  then  produce  a  sequence  of  points  along  that 
path.  In  general,  it  is  not  easy  to  develop  a  good  picture  of  a  multi¬ 
dimensional  manifold  from  information  along  one-dimensional  paths;  thus 
there  is  growing  interest  in  computational  methods  which  generate  multi¬ 
dimensional  grids  of  solution  points.  Up  to  now,  the  only  such  method 
appears  to  be  that  of  E.L.  Allgower  and  P.H.  Schmidt  [I].  It  utilizes  a 
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simplicial  continuation  algorithm  to  triangulate  a  p-dimensional  manifold 
Dy  means  of  p-simplices. 

in  [10]  a  new  algorithm  was  developed  for  computing  vertices  of  a 
triangulation  (Dy  p-simplices)  of  certain  suDsets  of  a  p-dimensional 
solution  manifold  (1.2).  It  depends  on  an  algorithm  for  constructing  a 
moving  frame  on  these  suDsets  of  M.  We  present  here  an  overview  of 
these  two  algorithms  and  illustrate  their  effectiveness  with  some 
numerical  examples. 


2.  Local  Coordinate  Systems 

At  any  point  x  of  M  the  tangent  space  TXN  may  be  identified  with 
the  kernel  of  the  Jacobian  DF(x), 

Tx n  =  ker  DF(x)  =  {  u  £  Rn  ;  DF(x)u  =  0  },  (2.  i ) 

and  then  the  corresponding  normal  space  Nxn  is  specified  as  tne 
orthogonal  complement  N*M  =  T/M-  =  r ge  DF(x)f. 

A  given  p-dimensional  subspace  T  c  Rn  induces  a  local  coordinate 
system  of  M  at  any  point  x  e  n  where 

T  n  NXI1  =  {0}  (2.2) 

As  shown,  for  instance,  in  [4]  or  [9],  at  any  x  e  M  where  (2.2)  holds  there 
exist  neighborhoods  V!  c  T  and  V2  £  Rn  of  the  origins  of  T  and  Rn. 
respectively,  and  a  unique  Cr'1  function  w:  Vj  -  T1,  w(0)  =  0.  such  that 

n  n  v2  =  {  y  £  Rn;  y  =  x  *  t  *  w(t),  t  e  Vj  }  (2.3) 

A  well-known  procedure  for  computing  tangent  bases  is  provided  by 
the  OR-decomposition 


DF(x)T  =  Q  [R]  .  Q  =  (Q i ,Q2)  . 

LoJ 


(2.-4) 


where  the  n  x  n  matrix  Q  is  orthogonal.  Q!  has  m  columns,  and  the  m  <  m 
matrix  R  is  upper  triangular  and  non-singular  for  x  £  n.  Then  the  p 
columns  of  Q2  form  an  orthonormal  basis  of  Txtt. 

If  x  e  n  is  a  point  where  the  QR- decomposition  (2.4)  has  been 
computed,  then  with  any  starting-point  y  =  y°  sufficiently  near  x  in  x*Tyn 
we  may  apply  the  chord-Gauss-Newton  process: 

For  k=0.l....  until  convergence 

1)  solve  RTz  =  F(y)  for  z  £  RP  (2.5) 

2)  y  ==  y  -  Q(z.0)T 

The  convergence  theory  of  these  methods  is  well  understood.  In  particular, 
a  theorem  of  Deuflhard  and  Heindl  (3)  can  be  used  to  ensure  that  there 
exists  for  any  x  £  n  a  neighborhood  V  =  V(x)  of  x  m  x-Txn  such  that  for 
any  y  in  v(x)  the  process  (2.4)  converges  to  some  y*  £  M.  Moreover,  we  can 
show  readily  that  y*-y°  e  NXM  and  hence  that,  in  the  notation  of  (2.3).  we 

have  y*  =  x*t+w(t),  t=y°-x.  In  other  words,  the  process  (2.5)  represents  an 
implementation  of  the  “corrector*  mapping  w  of  the  local  coordinate 
representation  (2.3). 

Thft  Moving  Frame  Algorithm 

Recall  that  a  vector  field  of  class  Cs,  s  s  r.  on  an  open  subset  M0  of 
M  is  a  Cs  function  u  :  M0  -*  TM  into  the  tangent  bundle  TM  such  that  u(x) 
belongs  to  TXM  for  x  each  M0.  A  moving  frame  of  class  Cs  on  M0 

associates  with  each  x  of  M0  an  ordered  basis  (frame)  (u1 . uPl  of  txm 

such  that  each  coordinate  map  u'  :  Mq  -»  TM,  1=1 . p  defines  a  vector  field 

of  class  Cs  on  M0.  We  shall  consider  only  orthonormal  moving  frames. 


In  our  setting,  an  algorithm  for  constructing  a  moving  frame  has  to 
generate  for  each  x  of  some  open  subset  n0  of  M  an  n  x  p  matrix  iKx) 
with  orthonormal  columns  such  that  DF(xXXx)  =  0  and  that  the  mapping 
U  :  M0  -*  RP*n  is  of  class  Cs  on  n0.  as  noted  in  [2],  the  QR-decomposition 
(2.4)  does  not  produce  continuously  varying  matrices  u(x).  This 
observation  extends  to  other  algorithms  of  a  similar  nature.  The  three 
remedies  proposed  in  [2]  do  not  concern  the  generation  of  a  moving  frame. 

For  the  moving  frame  algorithm  developed  in  [10]  we  assume  that 
some  method  is  available  for  computing  at  the  points  x  of  n  some  n  x  p 
matrix  Uq(x)  with  orthonormal  columns  that  span  Txn.  Of  course.  U0(x)  is 
not  expected  to  depend  continuously  on  x.  For  instance,  we  may  use  the 
QR-decomposition  (2.4). 

The  algorithm  is  based  on  the  selection  of  annxp  reference  matrix 
Tr  with  orthonormal  columns.  Then  for  a  point  x  of  the  manifold  we 
proceed  as  follows: 

(1)  Compute  the  tangent  basis  matrix  Uo(x) ; 

(2)  form  U0  :=  U0(x)Tr  ; 

(3)  compute  the  singular  value  decomposition  (3. 1 ) 

ATU0B  =  I  and  save  A  and  B  ; 

(4)  with  Q  =  ABT  form  the  basis  matrix  Uq(x)Q. 

The  following  result,  proved  in  [10],  guarantees  the  validity  of  this 
algorithm-- 

Theorem;  Let  M0  be  the  open  subset  of  M  where  the  subspace  of  Rn 
spanned  by  the  columns  of  the  reference  matrix  Tr  induces  a  local 
coordinate  system.  Then  the  mapping  x  z  M  *  Uo(x)Q  z  RnxP  given  by  the 

algorithm  (3.1)  is  of  class  Cr'1  on  Mo  and  defines  an  orthonormal  moving 
frame  on  M0. 

If  the  QR-decomposition  is  used  in  step  (!)  and  the  dimension  of  the 
manifold  is  small  in  comparison  with  the  space  dimension,  then  the 


principal  cost  of  (3.1)  derives  from  the  approximately  (2/3)n3  flops 
needed  for  the  decomposition  of  DF(x)T. 


In  practice,  it  has  turned  out  to  be  advantageous  to  construct  the 
reference  matrix  Tr  in  the  following  manner.  We  select  a  reference  point 

xr  on  n  Then  the  Euclidean  norms 

fj  =  |  Uo(xr)T  e*  |2  .  i=l . n 

of  the  rows  of  Uo(xr)  are  the  cosines  of  the  principal  angles  between  the 
tangent  space  of  M  at  xr  and  the  i-th  natural  basis  vector  e1  of  Rn.  The  r, 

are  independent  of  the  choice  of  the  basis  matrix  Uq(x).  Let  i[ ip  be 

the  indices  of  the  p  largest  of  these  tj  (with  ties  broken,  say. 
lexicographically).  Then  we  form  the  desired  reference  matrix  Tr  as  the 
matrix  with  the  columns  e1  .  i=ilt...,iD.  This  construction  is  analogous  to 
the  local  parameter  selection  in  the  continuation  program  PITCON,  [II]. 


±.  The  Tri angulation  Algorithm 

For  the  triangulation  of  a  p-dimensional  manifold  we  begin  by 
constructing  a  reference  triangulation  on  RP.  Let  I  be  the  collection  of 
simplices  of  this  triangulation.  Except  for  considerations  of 
computational  efficiency  and  simplicity,  no  restrictions  are  placed  on  z. 
We  refer,  for  example,  to  [12]  for  various  algorithms  for  triangulating  RP. 
For  our  purposes,  the  well-known  Kuhn-triangulations  have  been  useful. 

and.  in  the  case  p  =  2.  triangulations  of  R2  by  equilateral  triangles  have 
been  applied  as  well. 

Let  £  denote  a  given  vertex  of  this  triangulation  in  RP  and  h  >  0  a 
fixed  steplength.  Then  for  any  point  x  s  n  where  a  basis  matrix  u  of  TyM 
is  known,  the  mapping 


transfers  £  from  RP  onto  x  ♦  Txfl.  As  before,  let  V(x)  c  x  ♦  TXM  denote 
the  local  convergence  domain  of  the  Gauss-Newton  process  (2.5) .  If  ti  is 
a  vertex  of  Z  for  which  At\  e  V(x),  then  (2.5)  can  be  used  to  map  A-q  into 
a  point  yen.  The  set  TCS.x.u)  of  vertices  of  Z  that  can  be  mapped  onto  M 
in  this  way  shall  be  called  the  "patch”  corresponding  to  S.x.U  .  (The 
steplength  h  will  be  held  fixed  throughout). 

An  "idealized"  form  of  our  algorithm  can  now  be  phrased  as  follows: 

(1)  Select  a  reference  vertex  £*  of  z  : 

(2)  Select  a  reference  point  x*  e  M  and  let  Mq  be  the  subset 

where,  by  the  theorem,  the  moving  frame  algorithm  applies  ; 

(3)  Set  x  =  x*.  $  =  ; 

(4)  nark  the  vertex  £  as  “used" ; 

(5)  While  x 

(5  «  t  as  a  "center" 

(5jJ  Compute  the  frame  Ux)  by  algorithm  (3.1) ; 

(5c)  Select  all  vertices  of  the  patch  r(£,x,U(x)) 
which  have  not  yet  been  marked  "used”  ; 

(5d)  nap  these  vertices  onto  n  and  mark  them  “used”  ; 

(5e)  Choose  a  "used"  vertex  £  of  Z  not  marked  a  "center" 
and  let  x  be  its  computed  image  on  n  ; 

The  points  computed  on  n  inherit  the  connectivity  pattern  of  the 
original  simplices  of  Z  which,  in  turn,  induces  a  simplicial  approximation 
nr  of  n  in  Rn. 

The  algorithm  is  still  "idealized"  because,  in  practice,  it  is 
impossible  to  check  the  condition  x  e  n0  and  to  identify  the  vertices  of  z 
that  belong  to  r($,x,U(x)).  Thus,  special  provisions  have  to  be  added  in 
order  to  overcome  the  possible  failures  due  to  these  missing  checks.  We 
shall  not  go  into  details  here.  The  principal  approach  is  to  select  a 
"standardized"  patch  of  Z  which  is  used  in  step  (5c)  in  place  of  r(S.x,U(x)). 
Then,  in  step  (5d),  appropriate  alternatives  are  introduced  for  all  vertices 
where  a  failure  of  the  corrector  iteration  is  encountered. 


As  noted,  for  two-dimensional  manifolds  a  reference  triangulation 
of  equilateral  triangles  can  be  used.  Then  the  “standardized”  patch  is  the 
hatched,  star-shaped  region  in  the  center  of  Figure  l.  At  each  vertex,  the 
second  of  the  two  integers  is  a  counter  and  the  first  one  identifies  the 
"center”  £  that  is  used  in  mapping  that  vertex  onto  M.  Thus,  after  the 

reference  vertex  0,  the  nodes  7 . 12  become  centers  which  serve  to  map 

the  nodes  13 42  onto  M.  Then  the  process  continues  with  nodes 

17,18,19,23.24,28,29,33,34,38,39,42  as  centers.  This  is  no  longer  shown 
in  the  figure,  but,  in  practice,  we  always  continued  through  this  further 
stage.  It  results  in  a  total  of  1 14  triangles  on  M  and  involves  19  centers 
and  hence  as  many  Jacobian  evaluations.  This  indicates  the  efficiency  of 
the  algorithm.  In  fact,  in  terms  of  computed  points  per  Jacobian 
evaluation,  the  method  performes  better  than  most  continuation  processes. 
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Figure  1 


We  present  now  a  few  numerical  examples  to  indicate  the 
performance  of  the  methods.  But  space  limitations  force  us  to  be  brief 
flore  extensive  examples  will  be  given  elsewhere. 

Our  first  example  concerns  the  well-known  Belousov-Zhabotinskn 
reaction  [13].  As  in  [6]  we  write  the  mass  balance  equations  in  the  form 

(P-Xl)*2+  Xi(l-Xi)  -  EiflXi  =0 

-(M+x1)x2  ♦  x3  ♦  e20(o<-X2)  =0  (5.1) 

x  i  -  Xj(1  -$)  =0 

If  £j  =  1/1,500.  e2  =  1/56,250,  and  ji  =  8.4x1  O'6,  then,  as  discussed  in 
[6],  there  is  an  isola  point  approximately  at  the  point  with  the  coordinates 


Kl  =  0.249,  X2  =  0.750,  X3  =  0.125,  cx  =  3.508,  $>  =  0.997.  This  point  w3$ 
used  as  our  reference  point  on  n.  and  Figure  2  snows  tne  computed 
simp  1  icial  triangulation  (based  on  the  reference  triangulation  of  Figure  I). 
The  printed  page  is  the  cxj-plane  and  X2  is  the  third  coordinate  in  the 
figure. 

Our  second  example  concerns  the  roll  stability  of  maneuvering 
airplanes.  Without  going  into  details,  we  use  the  equations  originally 
formulated  in  [7]  and  given  in  [8]  and  [5]  in  a  simp!  if  led  form  Ax+$(x)  =  0, 
x  e  R8  .  Here  A  is  a  5  x  8  matrix  and  <t>:  R8  -*  R5  a  quadratic  function.  The 
(dimensionless)  control  parameters  x6,x7,x8  denote  the  elevator,  aileron, 
and  rudder  deflections,  respectively.  The  bifurcation  diagram  for  rudder 
deflections  x0  =  0  was  given  in  [8]  and  again  (with  some  extensions)  in  (5). 
In  the  neighborhood  of  the  origin  of  the  x6.x7-plane  it  has  the  form  shown 
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Figure  4 

The  examples  indicate  that  the  algorithms  work  very  efficiently, 
even  around  singularities.  Thus,  as  intended,  they  do  indeed  provide  a  new 
tool  for  deriving  information  about  the  shape  and  features  of  the  manifold. 
Of  course,  besides  any  graphical  representation,  the  extensive  numerical 


output  of  the  process  contains  a  wealth  of  further  information.  For 
instance,  linear  interpolation  between  the  computed  points  defines  the 
earlier  mentioned  simplicial  approximation  Mc  of  M.  The  corrector  process 
can  be  started  from  any  point  of  Ms  to  produce  additional  points  of  M.  In 
addition,  for  any  given  functional  it  is  easy  to  compute  a  contour  plot  of 
its  values  on  m£.  For  instance,  in  some  structural  problems  it  may  be  of 
interest  to  determine  lines  of  constant  stress  components.  Similarly,  the 
foldlines  on  II  represent  contour  lines  with  respect  to  a  measure  of  the 
orientation  of  the  projection  of  the  tangent  spaces  onto  the  parameter 
space.  This  provides  for  a  simple  method  of  approximating  the  fold-lines 
on  n  which  can  then  be  used  to  compute  the  fold  points  themselves  by 
means  of  one  of  the  numerous  local  iterative  processes  available  for  that 
purpose.  Examples  of  these,  and  other  post-processing  procedures  will  be 
given  elsewhere. 
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